setwd("/home1/30/jc227089/SRE_inv/ran_outs")

d<-c()
for (i in 1:50){
	fname<-paste("SRE_inv", i, ".RData", sep="")
	load(fname)
	d<-rbind(d, sim.result)
	rm(sim.result)	
}
save(d, file="simconcatDC01.RData")
setwd("/home1/30/jc227089/SRE_inv/Figures")
diff<-d[,3]-d[,2]
diff<-ifelse(is.na(diff), 3000, diff)  #treat NAs as 3000 for calculation of mean
d<-cbind(d, diff)
dsumm<-aggregate(d, by=list(d[,1]), FUN=mean)
dsumm.sd<-aggregate(d, by=list(d[,1]), FUN=sd)

X<-as.numeric(as.character(rownames(dsumm))) 
png("SRE_invFigDC01.png", width=10, height=20, units="cm", res=300)
par(cex.lab=1.5, mar=c(3,5,2,2), mfrow=c(3, 1), oma=c(3, 0, 0, 0))
sset<-cbind(dsumm$var.h1, dsumm$var.h2, dsumm$var.h3)
var.poly<-c(dsumm$var.h3-dsumm.sd$var.h3, rev(dsumm$var.h3+dsumm.sd$var.h3))
ymin<-min(sset, na.rm=T)
ymax<-max(c(max(var.poly, na.rm=T), max(sset, na.rm=T)))
matplot(X, as.matrix(sset), type="l", lty=1:3, col="black", xlab="", ylab=expression("Variance at " * italic("H")), bty="l", ylim=c(0, 1.05*ymax))
over.col<-col2rgb("grey70") #convert colour to rgb with 50% transparency
over.col<-rgb(red=over.col[1], green=over.col[1], blue=over.col[1], alpha=125, maxColorValue=255)
polygon(x=c(X, rev(X)), y=var.poly, col=over.col, border=NA)
legend(X[32], 0.9*ymax, legend=c("Initial", "Post-spread", "At breach"), lty=1:3, bty="n")
text(X[1], 1.05*ymax, labels="(a)")
sset<-cbind(dsumm$mean.d1, dsumm$mean.d2, dsumm$mean.d3)
var.poly<-c(dsumm$mean.d3-dsumm.sd$mean.d3, rev(dsumm$mean.d3+dsumm.sd$mean.d3))
ymin<-min(sset, na.rm=T)
ymax<-max(sset, na.rm=T)
matplot(X, as.matrix(sset), type="l", lty=1:3, col="black", xlab="", ylab="Mean dispersal probability", bty="l")
polygon(x=c(X, rev(X)), y=var.poly, col=over.col, border=NA)
legend(X[32], 0.9*ymax, legend=c("Initial", "Post-spread", "At breach"), lty=1:3, bty="n")
text(X[1], ymax, labels="(b)")
plot(dsumm$V1, dsumm$diff, bty="l", xlab="", ylab="Time to breach barrier (gens)")
lines(lowess(dsumm$V1, dsumm$diff))
text(min(dsumm$V1), max(dsumm$diff), labels="(c)")

#arrows(dsumm$V1, dsumm$diff-dsumm.sd$diff, dsumm$V1, dsumm$diff+dsumm.sd$diff, length=0)
#abline(reg=lm(dsumm$diff~dsumm$V1))
mtext("Distance of spread", side=1, line=1, outer=T)
dev.off()